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Abstract 

In this paper, we propose a new numerical method to compute the ground 
state solution of trapped interacting Bose-Einstein condensation (BEC) at 
zero or very low temperature by directly minimizing the energy functional 
via finite element approximation. As preparatory steps we begin with the 
3d Gross-Pitaevskii equation (GPE), scale it to get a three-parameter model 
and show how to reduce it to 2d and Id GPEs. The ground state solution is 
formulated by minimizing the energy functional under a constraint, which is 
discretized by the finite element method. The finite element approximation 
for Id, 2d with radial symmetry and 3d with spherical symmetry and cylindri- 
cal symmetry are presented in detail and approximate ground state solutions, 
which are used as initial guess in our practical numerical computation of the 
minimization problem, of the GPE in two extreme regimes: very weak inter- 
actions and strong repulsive interactions are provided. Numerical results in 
Id, 2d with radial symmetry and 3d with spherical symmetry and cylindrical 
symmetry for atoms ranging up to millions in the condensation are reported to 
demonstrate the novel numerical method. Furthermore, comparisons between 
the ground state solutions and their Thomas-Fermi approximations are also 
reported. Extension of the numerical method to compute the excited states 
of GPE is also presented. 
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1 Introduction 

Recently, there have been experiments of Bose-Einstein condensation (BEC) in di- 
lute bosonic atoms (alkali and hydrogen atoms) employing magnetic traps at ultra- 
low temperatures [7, 13, 23]. The condensation can consist of few thousand to mil- 
lions of atoms confined by the trap potential. This peculiar state of matter, whose 
existence was postulated back in the 1920s by Bose [12] and Einstein [22], exhibits 
several characteristics which set it apart from other homogeneous condensed matter 
systems [30, 18]. In fact, besides internal interactions, the macroscopic behavior 
of BEC matter is highly sensitive to the shape of the external trapping potential. 
Theoretical predictions of the properties of BEC matter can now be compared with 
experimental data by adjusting some tunable external parameters, such as the trap 
frequency and/or aspect ratio. Needless to say, these dramatic progresses on the 
experimental front are stimulation a corresponding wave of activity on both the 
theoretical and the numerical fronts. 

The properties of a BEC at temperatures T very much smaller than the critical 
temperature T c [30, 25] are usually described by the nonlinear Schrodinger equation 
(NLSE) for the macroscopic wave function known as the Gross-Pitaevskii equation 
(GPE) [27, 36]. Note that equations very similar to the GPE also appear in nonlinear 
optics where an index of refraction which depends on the light intensity leads to a 
nonlinear term like the one encountered in the GPE. 

There has been a series of recent studies which deal with the numerical solution 
of the time-independent GPE for ground state and the time-dependent GPE for 
finding the dynamics of a BEC. For numerical solution of time-dependent GPE, Bao 
et al. [9, 10, 11] presented a time-splitting spectral method, Ruprecht et al. [37] and 
Adhikari et al. [4, 5] used the Crank-Nicolson finite difference method to compute 
the ground state solution and dynamics of GPE, Cerimele et. al. [15] proposed a 
particle-inspired scheme. For ground state solution of GPE, Edwards et al. [21] 
presented a Runge-Kutta type method and used it to solve Id and 3d with spherical 
symmetry time-independent GPE. Adhikari [2, 3] used this approach to get the 
ground state solution of GPE in 2d with radial symmetry. Other approaches include 
an explicit imaginary-time algorithm used by Cerimele et al. [14] and Chiofalo et 
al. [16], a direct inversion in the iterated subspace (DIIS) used by Schneider et al. 
[38], and a simple analytical type method proposed by Dodd [20]. 

In this paper, we propose a new numerical method to compute the ground state 
solution of Bose-Einstein condensation by directly minimizing the energy functional 
through the finite element discretization. We begin with the 3d Gross-Pitaevskii 
equation, make it dimensionless to obtain a three-parameter model, show how to 
approximately reduce it to a 2d GPE and a Id GPE in certain limits. The ground 
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state solution is formulated by directly minimizing the energy functional under a 
constraint. Furthermore we present in detail the finite element approximation of 
the minimization problem for Id, 2d with radial symmetry and 3d with spherical 
symmetry and cylindrical symmetry (this is the most popular case in the current 
experiments of BEC), and provide approximate ground state solutions, which are 
used as initial guess in our practical numerical computation, in two extreme regimes: 
very weak interactions and strong repulsive interactions. Numerical results in Id, 2d 
with radial symmetry and 3d with spherical symmetry and cylindrical symmetry for 
atoms ranging up to millions in the condensation are reported to demonstrate the 
numerical method. Furthermore, comparisons between the ground state solutions 
and their Thomas-Fermi approximations are also reported by using our numeri- 
cal method. Our numerical results show that the Thomas-Fermi approximation 
is a 'good' approximation to the ground state solution in the strong repulsive in- 
teraction regime, but a 'worse' approximation in the medium interaction regime. 
Convergence rate of the Thomas-Fermi approximation to the ground state solution 
as a function of the number of atoms in the condensation is also reported by our 
method. Furthermore, we also extend our method to compute the excited states of 
GPE in Id. 

The paper is organized as follows. In Section 2 we begin with the 3d GPE, 
scale it to get a three-parameter model, show how to reduce it to lower dimensions. 
In Section 3 we present a new method to compute the ground state solution of 
a BEC by directly minimizing the energy functional and provide the approximate 
ground state solution in two extreme regimes: very weak interactions and strong 
repulsive interactions. In Section 4 we present detailed numerical formula and its 
finite element approximation for the ground state solution of GPE in Id, 2d with 
radial symmetry and 3d with spherical symmetry (in these cases, it is reduced to 
a Id problem), and in section 5 for 3d with cylindrical symmetry (in this case, it 
is reduced to a 2d problem). In Section 6 we report on numerical results of the 
ground state solution of BEC in Id, 2d with radial symmetry and 3d with spherical 
symmetry as well as cylindrical symmetry. In Section 7 some conclusions are drawn. 



2 Gross-Pitaevskii equation 

In this section, we will present the Gross-Pitaevskii equation in three dimension, 
how to scale it to a three-parameter model and reduction to lower dimension. 

At zero or very low temperature, a BEC is well described by the macroscopic 
wave function tp(x.,t) whose evolution is governed by a self-consistent, mean field 
nonlinear Schrodinger equation known as the Gross-Pitaevskii equation [27, 36, 25]. 
If a harmonic trap potential is included, the equation becomes 

ijfWzA = -|^v^(x,t)+^ (uy + uy + uy) ^x,o+nw(x,oivcm), 

(2-1) 
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where x = (x, y, z) T is the spatial coordinate vector, m is the atomic mass, h = 1.05 x 
10~ 34 [J s] is the Planck constant, N is the number of atoms in the condensation, and 
u x , ojy and uj z are the angular trap frequencies in x, y and ^-direction, respectively. 

For the following we assume (w.r.o.g.) oo x < u y < uo z . When UJ X UJy 

the trap potential is isotropic. Uq describes the interaction between atoms in the 
condensation and has the form 

U = (2.2) 
m 

where a is the s-wave scattering length (positive for repulsive interaction and neg- 
ative for attractive interaction). It is necessary to ensure that the wave function is 
properly normalized. Specifically, we require 



/ 



>(x,t)| 2 dx= 1. (2.3) 

A typical set of parameters used in current experiments with 87 Rb is 

m = 1.44 x 10~ 25 [kg], uj x = u y = u z = 20n [1/s], 
a = 5.1 \nm] = 5.1 x 1(T 9 [ml, N : 10 2 ~ 10 7 . 



2.1 Dimensionless GPE 

Following the physical literatures [21, 28, 38, 19, 24], in order to scale the equation 
(2.1) under the normalization (2.3), we introduce 



t = uj x t, x=— , 4>(x,t) = ao /2 V>(x,t), with a = J— — , (2.4) 
a V 

where ao is the length of the harmonic oscillator ground state. In fact, here we 
choose 1/uj x and a as the dimensionless time and length units, respectively. Plug- 
ging (2.4) into (2.1), multiplying by — ttt2, an d then removing all ~, we get the 

following dimensionless Gross-Pit aevskii equation under the normalization (2.3) in 
three dimension 

i = -^(x, t) + V(xMx, t) + K |V(x, t) | 2 ^(x, t), (2.5) 

where 

1 / 2 , 22, 2 2\ u y u z U N 4iraN 

2 v y ' oj x uj x a%n<jJ x a 

If we plug the typical set of parameters into above parameters, we find the values 

a w 0.3407 x 10~ 5 [m], k ss 0.01881A^ : 1.881 - 188100. 

There are two extreme regimes: one is when k = o(l), then equation (2.5) 
describes a weakly interacting condensation. The other one is when then (2.5) 

corresponds to a strongly interacting condensation or to the semiclassical regime. 
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2.2 Reduction to lower dimension 

In the following two cases, the 3d Gross-Pitaevskii equation (2.5) can approximately 
be reduced to 2d or even Id [29, 31]. In the case (disk-shaped condensation) 

Ld x W Uy, tu z > U) x , Jy ps 1, 7 Z > 1, (2.6) 

the 3d GPE (2.5) can be reduced to 2d GPE with x = (x, y) T by assuming that the 
time evolution does not cause excitations along the z-axis since they have a large 
energy of approximately hu z compared to excitations along the x and y-axis with 
energies of about hu x . Thus we may assume that the condensation wave function 
along the z-axis is always well described by the ground state wave function and set 
[29, 31] 

i(> = ifa(x,y,t)il> 3 (z) with ip 3 (z) = ^ 2 \<f)g(x,y,z)\ 2 dxdy^J , (2.7) 

where <p g (x,y,z) (see detail in (3.6)) is the ground state solution of the 3d GPE 
(2.5). Plugging (2.7) into (2.5), then multiplying by ip 3 {z) (where /* denotes the 
conjugate of a function /), integrating with respect to z over (— oo, oo), we get 



. <9^ 2 (x,t) 
where 



= -iv 2 ^ 2 + \ (x 2 + 7,V + C) ^ + (« 1>i(z) dz) \H 2 ^2, (2.8) 



/oo rc 
z 2 \Uz)\ 2 dz+ / 
-oo J — < 



<#3' 2 



dz 



dz. 



■ C't 

Since this GPE is time-transverse invariant, we can replace ip2 — > $ which 
drops the constant C in the trap potential and obtain the 2d GPE, i.e. 

i = -\ V 2 ^ + \ [x 2 + 7 2 y 2 ) 4, + (« jT ^) dz) (2.9) 

The observables are not affected by this. 

In the case (cigar-shaped condensation) [29, 31] 

uj y ^>uj x , u z ^>u x , <t=^ jy > 1, 7z > 1, (2-10) 

the 3d GPE (2.5) can be reduced to Id GPE with x = x. Similarly to the 2d case, 
we derive the Id GPE [29, 31] 

. di){x, t) = + + ^ K J^^ 3 (y : z) dydz) \tp(x,t)\ 2 ip(x,t), 

(2.11) 

where ^ 

My,z) = |0 9 (x,i/,z)| 2 dx) , (2.12) 
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with <fi g (x,y,z) (see (3.6)) the ground state solution of (2.5). 

In fact, the 3d GPE (2.5), 2d GPE (2.9) and Id GPE (2.11) can be written in a 
unified way [31] 

i^^ = ~VV(x,0 + V d (x)^(x,0 + «d|^(x,OlV(x,0, x eM d , (2.13) 



dt 



where 



( 1 



k / ipt(z) dz, 



V d (x) = < 



K, 



1 / 2 , 2 2 i 



2 2 N 

7** i , 



The normalization condition to (2.13) is 

/ M(x,t)| 2 dx= 1. 



d= 1, 
d = 2, 



d = 3. 



(2.14) 
(2.15) 



3 Ground state solution 

In this section, we will propose a new numerical method by directly minimizing 
the energy functional via finite element discretization to compute the ground state 
solution of a BEC (2.13). Furthermore we will also provide approximate ground 
state solutions in two extreme regimes: very weak interactions and strong repulsive 
interactions. 



3.1 Minimization problem 

To find a stationary solution of (2.13), we write 

#M) = e-«"*#x), (3.1) 

where fi is the chemical potential of the condensation and <fi a real- valued function 
independent of time. Inserting into (2.13) gives the following equation for 0(x) 

H 0(x) = -iv 2 0(x) + V d (x)0(x) + K d |0(x)| 2 0(x), x G R d , (3.2) 

under the normalization condition 

/ |0(x)| 2 rfx=l. (3.3) 
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This is a nonlinear eigenvalue problem under a constraint and any eigenvalue \i can 
be computed from its corresponding eigenfunction by 



Jr< 



|V</.(x)| 2 + ^(x)0 2 (x) + ^ 4 (x) 



dx. 



(3.4) 



The Bose-Einstein condensation ground-state wave function </> s (x) is found by solv- 
ing this eigenvalue problem under the normalization condition (3.3) with the mini- 
mized chemical potential \x g . Usually, the ground state problem is formulated vari- 
ationally. Define the energy functional 



EM := f 



l -w<p\ 2 +v d {x) i0i 2 +|i0r 



d: 



'x. 



(3.5) 



It is easy to see that critical points of the energy functional E K (<f>) under the con- 
straint (3.3) are eigenf unctions of the nonlinear eigenvalue problem (3.2) under the 
constraint (3.3) and versus versa. In fact, (3.2) can be viewed as the Euler-Lagrange 
equation of the energy functional E K (<f>) under the constraint (3.3). To compute the 
ground state <f> g , we solve the minimization problem 

(V) Find (fig, 4> g G V) such that 

E K (<j) g ) = min E K {<j>\ \i g = fi K ((f) g ) = E K (<j> g ) + / — 0j(x) rfx, (3.6) 

where the set V is defined as 
V 



— |^ I E K (4>) < oo, j Rd |0(x)| 2 rfx=l|. 



In non-rotating BEC, the minimization problem (3.6) has a unique real valued 
nonnegative ground state solution g (x) > for x e M. d [32]. In physical literatures, 
the minimizer of (3.6) was obtained by either the Runge-Kutta type method [21, 2, 3] 
or the imaginary time method [6, 16, 14]. Here we present a method by directly 
minimizing the energy functional E K (<f)) through the finite element discretization 
[17]. 



3.2 Approximation in a bounded domain 

The eigenvalue problem (3.2) and the minimization problem (3.6) are defined in 
M. d . In practical computation, usually they are approximated by problems defined 
on a bounded computational domain. Since the full wave function must vanish 
exponentially fast as |x| — > oo and due to symmetry, choosing Ri,---,Ra > 
sufficiently large and denoting 

Q R = [0,^1] x ••• x [0,R d ], 

then the minimization problem (3.6) can be approximated by 
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/n« 2 



(V R ) Find (/if, 0f G K,) such that 

^W) = min £*(0), ^ = £*(«/>£) + jf 
where the set V g and the functional E^(4>) are defined as 

^ 5 = {</> | £*(</>) < oo, 2 d J^ R |0(x)| 2 dx=l, 0(i2i,x 2 ,...,x d ) 



(3.7) 



^(Xi,---,^-!,^) = 0>, 



E«{4>) 



-|V0(x)| 2 + ^(x)|0(x)| 2 + ^ |0(x)| 4 



dx. 



3.3 Discretization 

The functional in (3.7) (or £ K (0) in (3.6)) can be discretized by the finite 

element method [17], finite difference method [34] or spectral method [26]. Here 
we use the finite element method because it can easily keep the good properties of 
Eff(<f>), e.g. positive, coercive and weakly lower semicontinuous when K d > 0, on the 
unit sphere in finite dimensions. Let 

V g = {(f)\ E*(4>) < oo, (f)(R 1 ,x 2 , ■ ■ ■ , x d ) = ■ ■ ■ = 0(xi, • • • , x d _i, R d ) = O} 

and Vg be a finite dimensional subspace of V g [17], i.e. V g C V g . Then the finite 
dimensional set 



V 



e V! 



n R 



|0 h (x)| 2 dx 



is a subset of V g , i.e. V g C V g . Thus the minimization problem (3.7) can be 
approximated by 



{V R > h ) Find (nf h ,4>f h G V g h ) such that 

E«{<t>f h ) = min E«(d>% ^ = E%fi*) + L % 



t> h ev g h 



r K d 


>f fc (< 


Iqr 2 





dx. 



(3.8) 



Introducing a functional with a Lagrange multiplier corresponding to the nor- 
malization condition (3.3), i.e. 

E«{<t>\\) = ^^(^/^ |/(x)| 2 dx-l 
= F($ h ,A), <f) h eV g h , AgM, 

M 



where fe (x) = * h (x) • = 51 *i( x )$? with ^ ft (x) = (tff(x), • • • , *m( x )) a 
basis of the finite element subspace V g and $ h = ($i , • • ■ , the unknowns of 



(x) [17], then the minimizer 



LR,h 



^ (x) - §R' h of the minimization problem (3.8) 



is a critical point of the functional E^((f) , A). This implies 



V$ h F($ fc ,A) 



o, 



dF($ h ,\) 



= 0. 



(3.9) 



The nonlinear system (3.9) is solved by the Newton's method [35] or quasi-Newton's 
method [35] with a proper choice of the initial data (<pf h )^ = ^ A (x) • ($^)(°) and 
A(°) is the least square solution of 



V* h F($\A) 



^h=($«.h)(0) j A=A(0)) 



0. 



The initial guess ((f)^ h )^ is chosen as the interpoltant on V g h of the approximate 
ground state solution for very weak interactions (4.13) or strong repulsive interac- 
tions (4.15) when Kd is not too big or not too small, respectively. These approximate 
ground state solutions are given in the next subsection. Another way to choose the 
initial guess is to use a continuation technique, i.e. use the numerical solution of the 
ground state function for a small Kd as initial guess for computing the solution of a 
larger K d . 



3.4 Approximate ground state solution 

Here we present the approximate ground state solution of (3.2) in two extreme 
regimes: very weak interactions and strong repulsive interactions. These approxi- 
mate ground state solutions are used as initial guess (0^ ,ft )(°) in our practical com- 
putation of the minimization problem (3.8) ( or (3.9)). 

For a very weakly interacting condensation, i.e. Kd = o(l), we drop the nonlinear 
term (i.e. the last term on the right hand side of (3.2)) and get [31] 

H 0(x) = -^V 2 0(x) + V d (x) 0(x), x G R d . (3.10) 
The ground state solution of (3.10) is 




1 f e- 2 / 2 , d=l, 

C(x) 

I (7 w 7*) 1 / 4 e-( xa+ ™ a+ T , ** a >/ 2 , d = 3. 



(3.11) 

This can be viewed as an approximate ground state solution of (2.13) in the case of 
a very weakly interacting condensation. 

For strong repulsive interactions, i.e. Kd ^> 1, we drop the diffusion term (i.e. 
the first term on the right hand side of (3.2)) corresponding to the Thomas-Fermi 
approximation [31]: 

li 0(x) = \/ d (x)0(x) + K d |0(x)| 2 0(x), x G R d , (3.12) 
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The ground state solution of (3.12) is the compactly supported function <j) d (x): 

(^) 1/2 , d = 2, (3.13) 



[ 0, otherwise. 

Remark 3.1 As indicated in Figure 6, an interface layer correction has to be con- 
structed in order to improve the approximation quality in the Thomas-Fermi regime 
(i.e. Kd ^> I). For a convergence proof o/0| — > S as /c 3 — > +oo (without conver- 
gence rate) we refer to [32]. In Section 6, these convergence rates are reported based 
on our numerical solutions. 



4 Ground state solution in Id, 2d with radial sym- 
metry and 3d with spherical symmetry 

In this section, we present detailed numerical formula and its finite element approx- 
imation for the ground state solution of GPE in Id, 2d with radial symmetry (i.e. 
7 y = 1) and 3d with spherical symmetry (i.e. 7 y = j z = 1). In these cases, the 
problem is reduced to Id. Due to symmetry, the GPE (2.13) essentially collapses to 
a Id problem with r = |x| G [0, oo) for d — 1, 2, 3 

S = 0, ip(r, t) -> 0, when r -> oo. (4.2) 

or 

The normalization condition collapses to 

roo 

C d \1>(r,t)\ 2 r d - 1 dr=l, (4.3) 



o 

where 

'2, d=l, 
C d = I 2tt, d = 2, 
4tt, d = 3. 

4.1 Minimization problem 

The eigenvalue problem (3.2) collapses to 

^ <^ r ) = ~\^^ ( r "~ 1 ^ (r) ) + Y 0(r) + Kd W r )l^( r )' < r < oo, (4.4) 
0'(O) = 0, <j>(r) -> 0, when r -> oo, (4.5) 
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under the normalization condition 

C d / |0(r)| 2 r d ~ l dr = 1. 
Jo 

The minimization problem (3.6) collapses to 
(V) Find (fi g , <p g G V) such that 

E K (<P 9 ) = mm E K (<P), N = E K (<P g ) + C d f°° ^ <f> 4 g (r) r^ 1 dr, 

tp£V JO Z 

where the set V and the energy functional E K ((f>) are defined as 

V={</>\ EM < oo, C d jH |0(r)| 2 r d -' dr = l} , 



(4.6) 



(4.7) 



E K {4>) = C d 



„<2-l 



'(r)] 2 + r 2 W + «d 0» 



dr. 



4.2 Approximation in a bounded domain 

The eigenvalue problem (4.4), (4.5) and the minimization problem (4.7) are defined 
in a semi-infinite interval (0, oo). In practical computation, usually they are ap- 
proximated by problems defined on a finite interval. Since the full wave function 
must vanish exponentially fast as r — > oo, choosing R > sufficiently large, then 
the eigenvalue problem (4.4), (4.5) can be approximated by 



H 0(r) = - 



1 1 d 

2 r d-i dr 



_i d 
dr 



<t>{r) + 770(0 + K d |0(r)| 2 0(r), < r < R, (4.8) 



0'(O) = 0, <f>(R) = 0, 



under the normalization condition 

C d f R |0(r)| 2 r*- 1 dr = 1. 

Similarly the minimization problem (4.7) can be approximated by 
(V R ) Find (/if, 0f G V g ) such that 



(4.9) 



(4.10) 



R 



«) = min^), tf = E?(<l>?)+C d r% <f>f(r)V r d ~ l dr, (4.11) 



,l4 



o 2 



where the set V g and the functional E^((j>) are defined as 



= I £f (0) < 00, C d J q 

e«{4>) = a 



r r )|2 r d-i 



dr = 1, 0(i?) = L 



i? y.d—1 



'(r)] 2 + r 2 2 (r) + /t d </> 4 (r) 



dr. 
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4.3 Finite element approximation 

Assume that 

= r < r l < r 2 < ■ ■ ■ < r M = R 
is a partition of the interval [0, R] with mesh size h [17]. Let 

V g h = {4>\r) G C([0,R}) | h (r)| [rj , rj+l] G P 1 ([r j , r, +1 ]), < j < M - 1, 
/OR) = 0}, 

J? = j^^G^IC,^ |/(r)| 2 r^ 1 rfr = l,|; 

where Pi denotes piecewise linear polynomials. Then the finite element approxima- 
tion of the problem (4.7) is 

(V R ' h ) Find {nf h ,<f>f h e O such that 



R(±R,h\ 



KM 



R 



min^(0), p* h = E?(<ft*)+C d r% 



o 2 



bf h (r) r d - x dr. (4.12) 



4.4 Approximate ground state solution 

In these cases, the approximate ground state solution collapses to the following. For 
a very weakly interacting condensation, i.e. n d = o(l), the ground state solution is 



d 



CM = ie- r2 / 2 , d=l,2,3. 



2' rrfW 7T d / 4 

For strong repulsive interactions, i.e. Kd ^> 1, the ground state solution is 



((d+l) 2 -l)K d 

a. 



2/(d+2) 



d= 1,2,3, 



(4.13) 



(4.14) 



0, 



(^ g -ry2)/K d , r 2 <2^ d=1>2>3 . (4 . 15) 



otherwise, 



5 Ground state solution in 3d with cylindrical 
symmetry 

In this section, we present detailed numerical formula and its finite element approx- 
imation for the ground state solution of GPE in 3d with cylindrical symmetry (i.e. 
7 y = 1). In this case, the problem is reduced to 2d. Due to symmetry, the GPE 
(2.13) essentially collapses to a 2d problem with r = ^/x 2 + y 2 G [0, oo) 



djj(r,z,t) 
dt 



1 d ( dtp\ d 2 ip 
r dr \ dr j dz 2 



+ \(r 2 + 1 2 z z 2 )^ + K\^, (5.1) 
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dip(0,z,t) 



0, < z < oo, 



< r, z < oo, 

&i/>(r,0,t) 



dr dz 
ifj(r, z, t) — > 0, when r + |^| — > oo. 



0, < r < oo, 



The normalization condition collapses to 

/•OO /'OO 

47r / / Wr, z,t)\ 2 r drdz = 1. 
JO Jo 

5.1 Minimization problem 

The eigenvalue problem (3.2) collapses to 



H <f>(r, z) 
d<j>(0,z) 



ld_ ( df\ dV 

r dr \ dr J dz 2 

< r, z < oo 
d(j)(r, 0) 



+ i(r 2 + 7,V)0 + K |0p 



0, < z < oo, 



<9r cte 
<f>(r, z) — > 0, when r + |^;| — > oo, 



0, < r < oo, 



under the normalization condition 



47T 



JO 



(r, z)\ 2 r drdz = 1. 



The minimization problem (3.6) collapses to 
(V) Find (/i 3 , 4> g G V) such that 



E K ((f) g ) = min £ K (0), = E K (<p g ) + 4tt 



K 



o Jo 2 T9 



where the set V and the functional E K {(f>) are defined as 



V : 

EM = 47T 



| E K ((f>) < OO, 47T 

oo roo rp 

o 7o 2 



oo /-oo 



(5.2) 
(5.3) 



(5.4) 



(5.5) 

(5.6) 
(5.7) 

(5.8) 



- <t> q {r,z)r drdz, (5.9) 



|(/>(r, z)| r drdz = 

o Jo 

<f) 2 r {r, z) + 4> 2 z {r, z) + (r 2 + -fiz 2 )</> 2 (r, z) + k 4 (r, z) 



drdz. 



5.2 Approximation in a bounded domain 

The eigenvalue problem (5.5)-(5.7) and the minimization problem (5.9) are defined 
in the first quadrant of the rz-plane. In practical computation, usually they are 
approximated by problems defined on a bounded domain. Since the full wave func- 
tion must vanish exponentially fast as r + \z\ — > oo, choosing R > and Z > 
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sufficiently large, then the eigenvalue problem (5.5)-(5.7) can be approximated by 



H <f>(r, z) 

d_mz) 

dr 



1 d 



d 2 



dr \ dr I dz 2 



+ -(r 2 + 1 2 z 2 )<j) + K\<j)\ 2 <j>, 



= 0, < z < Z, 



< r <R, 0<z < Z, 
d<j>(r, 0) 



dz 



= 0, < r < R, 



4>{R,z) = Q, 0<z<Z, 
under the normalization condition 

cR rZ 



(r,Z)=0, < r < R, 



(5.10) 

(5.11) 
(5.12) 



4tt / / |0(r,z)| 2 r dzdr = 1. (5.13) 
Jo Jo 

Similarly the minimization problem (5.9) can be approximated by 
(V R ) Find (jif, (f)f G V g ) such that 



E^f) = min E*{4>), = E?{<%) + 4vr 



R rZ 



4>&v g " ' Jo Jo 2 

where the set V g and the functional J R {4>) are defined as 



(f>q(r,z)] rdrdz, (5.14) 



V g = {</> | £*(0) < oo, 4tt jf* |0(r,z)| 2 r dzdr = 1, 
0(i?, z) = 0, < z < Z, 0(r, Z) = 0, < r < 



'o jo 2 



r, z) + (j) z (r, z) + (r + 7 2 z )0 (r, z) + k (f) (r,z 



dzdr. 



5.3 Finite element approximation 

Assume that 

= r < r\ < r 2 < ■ ■ ■ < r M = R, = z < z\ < z 2 < ■ ■ ■ < Z N = Z 
is a partition of the rectangle [0, R] x [0, Z] with mesh size h [17]. Let 



yh 
9 



<l> h (r,z)eC([0,R}x[0,Z})\<f> h (r , Z)\[ rjt r j+1 ]x[zi,zi +1 ] e Q 1 ([r j ,r j+1 ] x [zi,zi +1 ]), 
< j < M - 1, < I < N - 1, <f) h (R, z) = 0, 0<z<Z, 
(j) h (r,Z) = 0, < r < i?J, 



V? 



G V7 | 4vr 



o Jo 



\<p h (r, z)\ 2 r dzdr = 1 
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where Qi denotes all bilinear polynomials. Then the finite element approximation 
of the problem (5.14) is 

(V R ' h ) Find (fif h , <pf h e V g h ) such that 

E^ h ) = mm £*(0), f,f h = E^f h ) + 4tt j* j* \ [(fif' h (r, z)}' r drdz. 

(5.15) 



5.4 Approximate ground state solution 

In this case, the approximate ground state solution collapses to the following. For a 
very weakly interacting condensation, i.e. k 3 = k = o(l), the ground state solution 
is 

1/4 

^ = 1 + y, <%(r,z) = ^j- 4 e- {r2+ ^ 2)/2 . (5.16) 
For strong repulsive interactions, i.e. Ks — k >> 1, the ground state solution is 
s 1/15K7A 275 

^ = 2(-srJ ' (5 - 17) 

#(r, z) = I V / ^-(- 2 + ^ 2 )/2)/-, - 2 + 7 Z V < 2/4, (51g) 
( 0, otherwise. 



6 Numerical results 

In this section we shall report on numerical error analysis of our method, numerical 
ground state solutions of (2.13) in Id, 2d with radial symmetry and 3d with spher- 
ical symmetry as well as cylindrical symmetry. Furthermore we also compare the 
numerical ground state solution of (2.13) in 3d and the corresponding Thomas-Fermi 
approximation (5.18). Convergence rates of the Thomas- Fermi approximations to 
their exact counterparts are also reported. 



6.1 Numerical error analysis 



In this subsection, we study numerically the convergence rate of the finite element 
discretization to the minimization problem (3.7) in Id. We choose d — 1 and 
&d = k>i = 62.742 in (4.11). We compute a numerical solution of (4.11) in Id on 
Q = [0,8] by using the discretization (4.12) with a very fine mesh, e.g. h = as 
the "exact" ground state solution 



Table 1 shows the errors E K ((fi g ) — E K {(fi g ), 



fig-fig, max 



X)-<j>g 



L 2 (ny 



- w 



L!(f2) 



and 



for different mesh size h. Here we use the standard Sobolev space norms [1]. 

From Table 1, we observe that the approximate energy E K ((f>g), chemical potential 
fig, ground state solution (fig, and atom density function ((fig) 2 converge to E K (<f> g ), 
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EM)-E K (</> 9 ) 



max 



S (x) 



/i 



L 2 (n) 
) 2 



6.528E-4 
3.462E-4 

3.222E-3 

2.177E-3 

1.042E-3 

2.493E-2 



1.570E-4 
8.073E-5 
8.450E-4 
5.290E-4 

2.633E-4 

1.246E-2 



3.878E-5 
1.986E-5 
2.091E-4 
1.323E-4 

6.579E-5 

6.218E-3 



9.560E-6 
4.887E-6 
5.300E-5 
3.394E-5 

1.700E-5 

3.091E-3 



2.274E-6 
1.160E-6 
1.359E-5 
9.295E-6 

4.927E-6 

1.508E-3 



Table 1. Numerical error analysis of the finite element discretization (4.12). 

/i g , 4> g in maximum norm and L 2 -norm, and (4> g ) 2 in Z^-norm, respectively, at second 
order convergence rate when the mesh size h goes to zero. Furthermore (f> g converges 
to (pg in i7 1 -norm at first order convergence rate. 

6.2 Results in Id, 2d with radial symmetry and 3d with 
spherical symmetry 

An interesting property of the condensation wave function in these cases is its root 
mean square (rms) size r rms defined by 

rLs = (r 2 )=C d r 2 <t> 2 {r)r d - l dr, d = 1,2,3. (6.1) 
Jo 

In our computations, we choose R = 16 in (4.8) with a uniform partition of the 
interval [0, R] of mesh size h = ^ in (4.12). 

Figure 1 shows the ground-state condensation wave function 4> g (r) (with (f> g (—r) = 
4> g (r) for r 6 R in Id) versus r and the chemical potential \x g for different Kd, and 
Table 2 lists 9 (O), r rms , /i g versus Kd for d — 1. For comparison, we also list the 
chemical potential \i g obtained by the Thomas- Femi approximation (TFA) in (4.14). 
Furthermore Figure 2 and Table 3 show similar results for d — 2, and Figures 3 and 
Table 4 for d = 3. 

From Figures 1-3 and Tables 2-4, we can see that the chemical potential fj, g and 
the root mean square size will increase when the interaction coefficient K d ( i.e. the 
number of atoms in the condensation) is increasing. On the other hand, the peak of 
the ground state solution (f) g (0) will decrease. If we use the typical set of parameter 
values in Section 2, a k 3 = 31371 corresponds to a condensation population of 
N 1 667 800, i.e. approximately one and a half millions atoms in the condensation. 
Furthermore the Thomas-Femi approximation gives accurate chemical potential \x g 
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0o(O) 






Numerical 


TFA 


-12.5484 


1.7718 


0. 


1444 


-19.669 


NA 


-6.2742 


1.2654 


0. 


2810 


-4.9553 


NA 


-2.5097 


0.9132 


0. 


5133 


-0.8061 


NA 





0.7511 


0. 


.7071 


0.5000 


NA 


3.1371 


0.6459 


0. 


.8960 


1.5265 


NA 


12.5484 


0.5297 


1. 


.2454 


3.5965 


3.538 


31.371 


0.4556 


1. 


.6416 


6.5526 


6.517 


62.742 


0.4060 


2. 


0495 


10.369 


10.345 


156.855 


0.3485 


2. 


.7679 


19.0704 


19.056 


313.71 


0.3105 


3. 


4823 


30.259 


30.249 


627.42 


0.2766 


4. 


.3847 


48.024 


48.018 


1254.8 


0.2464 


5. 


.5228 


76.226 


76.222 



Table 2. Ground state chemical potential n g , maximum value of the wave function 
4> g (0) and root mean square size r rms versus the interaction coefficient Kd in Id 
(d=l). 



«2 


0,(0) 


^rms 


Numerical 


TFA 


-5.8 


2.1770 


0.3208 


-5.552 


NA 


-4.5 


0.8948 


0.7091 


-0.2923 


NA 


-2.5097 


0.6754 


0.8775 


0.4997 


NA 





0.5642 


1.0000 


1.0000 


NA 


3.1371 


0.4913 


1.1051 


1.4200 


NA 


12.5484 


0.3919 


1.3068 


2.2558 


1.9986 


62.742 


0.2676 


1.7881 


4.6098 


4.4689 


313.71 


0.1787 


2.6044 


10.068 


9.9928 


627.42 


0.1502 


3.0845 


14.1892 


14.132 


1254.8 


0.1262 


3.6598 


20.0286 


19.9854 



Table 3. Ground state chemical potential [i g , maximum value of the wave function 
(f) g (0) and root mean square size r rms versus the interaction coefficient K d in 2d with 
radial symmetry (d — 2). 
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6 (0) 


T 

' rms 












IV n m pti or 1 


TFA 


7 


n 7R1 3 

U . 1 uio 






NA 


3 1 371 


U .4:00 1 


1 1 ^91 




NA 


n 

u 




1 99/L& 


1 ^nnn 

1 .JUUU 


N A 

IN 1 V 


3.1371 


0.3843 


1.2785 


1.6774 


NA 


12.5484 


0.3180 


1.3921 


2.0650 


1.4762 


31.371 


0.2581 


1.5356 


2.5861 


2.1298 


125.484 


0.1738 


1.8821 


4.0141 


3.7082 


627.4 


0.1066 


2.5057 


7.2484 


7.059 


3137.1 


0.0655 


3.4145 


13.553 


13.438 


31371 


0.0328 


5.3852 


33.810 


33.755 



Table 4. Ground state chemical potential maximum value of the wave function 
9 (O) and root mean square size r rms versus the interaction coefficient K d in 3d with 
spherical symmetry (d — 3). 



only when the interaction is very big and poor approximation for intermediate 
values of n d (cf. Tables 2-4). 

Here we also report the convergence rates of the Thomas-Femi approximations 
(p s d and fi s d to the exact ground state solution <p g and /i g , respectively, as ^ — > 0. 
Tables 5,6&7 list these results for d — 1, 2&3, respectively. 

From Tables 5-7, we observed numerically the convergence rates of the Thomas- 



Femi approximations as following: (1) in Id: \fi g — ~ 0(l/4' 57 )> 2 — (0f^ 2 
0(l/< 72 ), - 4>{\\ L2 « 0(l/< 43 ) and ||0 9 - 0-|| L4 » 0(l/< 41 ) as £ _> ; (2) in 

2d: \im 9 - i4\ * o(i/4 40 ), ||$ - (0!) 2 || L1 « o(i/4 58 ), - 0^ll L2 « o(i/4 31 ) 
and ||0 9 -^|| L 4 « 0(l/4 35 ) as ^ -> 0; and (3) in 3d: \fi g - fi 8 3 \ « 0(l/4 31 ), 
0(l/4 45 ), - 01|| L2 » 0(l/4 24 ) and ||0 9 - 0^|| L4 » 0(l/4 32 ) 



2 - 

as -> 0. 

K3 



L 1 



6.3 Results in 3d with cylindrical symmetry 

The interesting properties of the condensation wave function in this case are its root 
mean square (rms) sizes in r- and ^-direction r rms and z Tms , respectively, defined by 



r 2 

rms 



z 2 

rms 



poo poo 

= ( r 2 ) = 4vr / r 2 ^ 2 {r) r drdz, (6.2) 
Jo Jo 

poo roo 

= (z 2 ) = 4tt / / z 2 ip 2 {r) r drdz. (6.3) 

J0 J 



We present computations for two cases: 
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~ Mil 



1 

Ki 

Error 
Rate 



100 

1.875E-2 
0.5655 



200 

1.267E-2 
0.5711 



1 

400 

3.528E-3 
0.5775 



800 

5.715E-3 
0.5827 



1600 



3.816E-3 



<t>l - (0?) 2 



Error 9.915E-3 5.799E-3 3.451E-3 2.102E-3 1.319E-3 
l 1 Rate 0.7738 0.7488 0.7151 0.6723 



Error 5.709E-2 



Rate 



0.4346 



4.224E-2 3.120E-2 2.305E-2 
0.4371 0.4368 0.4375 



1.702E-2 



M B -<K\ 



L 4 



Error 
Rate 



5.676E-2 
0.4079 



4.278E-2 
0.4103 



3.219E-2 
0.4080 



2.426E-2 
0.4083 



1.828E-2 



Table 5. Convergence rates of the Thomas-Femi approximations as — — > in Id 

«i 

(d=l). 



1 

Error 
Rate 



200 

9.004E-2 
0.3947 



400 

6.849E-2 
0.4007 



1 

800 

5.188E-2 
0.4062 



1600 

3.915E-2 
0.4112 



3200 



2.944E-2 



Error 
Rate 



4.458E-2 
0.5835 



2.975E-2 
0.5881 



1.979E-2 
0.5941 



1.311E-2 
0.5828 



8.753E-3 



II || Error 1.292E-1 1.038E-1 

1109 Ml2 Rate 0.3158 0.3164 



8.336E-2 6.686E-2 
0.3182 0.3211 



5.352E-2 



Error 6.613E-2 5.169E-2 4.036E-2 
Rate 0.3554 0.3570 0.3580 



3.149E-2 2.449E-2 
0.3627 



Table 6. Convergence rates of the Thomas-Femi approximations as > in 2d 

(d = 2). 
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111111 









«3 


400 


800 


1600 


3200 


6400 


\fig - 






Error 
Rate 


2.169E-1 
0.3023 


1.759E-1 
0.3058 


1.423E-1 
0.3111 


1.147E-1 
0.3142 


9.225E-2 


«"( 


03) 2 




Error 
Rate 


1.054E-1 
0.4470 


7.732E-2 
0.4541 


5.644E-2 
0.4579 


4.109E-2 
0.4649 


2.977E-2 


11^- 


03ll L 2 


Error 
Rate 


2.046E-1 
0.2437 


1.728E-1 
0.2461 


1.457E-1 
0.2455 


1.229E-1 
0.2492 


1.034E-1 


u 9 - 




Error 
Rate 


6.524E-2 
0.3195 


5.228E-2 
0.3241 


4.176E-2 
0.3240 


3.336E-2 
0.3305 


2.653E-2 



Table 7. Convergence rates of the Thomas-Femi approximations as — = > in 

3d (d = 3). 

Case I. 87 Rb used in JILA with uj x —uj y < uj z [7]. The detailed data are 
m = 1.44 x 10~ 25 [kg], u x = u y = 20ir[l/s], u z = 4w x [l/s], 
a = J = 0.3407 x 10 -5 H, \a\ = 5A[nm], k = AnaN/a^ = 0.01881JV. 

V u x m 

Case II. 23 Na used in MIT (group of Ketterle) with oo x = oo y » uo z [8]. The detailed 
data are 

m = 3.8 x 10 _26 [A#], u x = u y = 720?r[l/s], u z = 7n[l/s], 

a = J— = 1.1209 x 10 _5 [m], \a\ = 2.75[wn], re = 47raN/a = 0.003083iV. 
V uj z m 



In case II, we choose a = ( m stead of a = \J ^-^) as cj z <C such 

that the root mean square size is of 0(1). The other parameters should be adjusted 
accordingly. 

Figure 4 shows the ground-state condensate wave function along r- and z-axis, 
4> g (r,0) and (f> g (0,z), respectively, for different re and surface plots of <f) g (r,z) for 
re = 15408 and re = 188.1, and Table 8 lists /i g , (f> g (0, 0), r rms , z rms versus k^ — k for 
case I. Figure 5 and Table 9 show similar results for case II. Furthermore Figure 6 
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a 


N 


K 


0,(0,0) 


^rms 


^rms 


Numerical 


TFA 


-5.1 [nm] 


250 


-4.705 


0.9926 


0.7468 


0.3268 


1.9294 


NA 


r 1 r -| 

-5.l|nmJ 


100 


-1.881 


0.6788 


0.9324 


0.3476 


2.7283 


NA 


5.l[nm] 








0.602 


1.000 


0.3539 


3.000 


NA 


5.l[nm] 


l 000 


18.81 


0.3824 


1.325 


0.3807 


4.362 


3.022 


5.l[nm] 


5 000 


94.05 


0.2477 


1.7742 


0.4214 


6.680 


5.752 


5.l[nm] 


10 000 


188.1 


0.2023 


2.041 


0.4497 


8.367 


7.591 


5.l[nm] 


50 000 


940.5 


0.1248 


2.842 


0.5532 


14.95 


14.45 


5.l[nm] 


100 000 


1881 


0.1012 


3.276 


0.6174 


19.47 


19.06 


5.l[nm] 


400 000 


7524 


0.0666 


4.341 


0.7881 


33.47 


33.19 


5.l[nm] 


800 000 


15048 


0.0540 


4.992 


8976 


44.02 


43.80 



Table 8. Ground state chemical potential fi g , maximum value of the wave function 
4> g (0, 0) and root mean square sizes r rms , z rras versus the interaction coefficient k 3 — k 
in 3d with cylindrical symmetry under case I. 



a 


N 


K 


0,(0,0) 


^rms 


^rms 


Numerical 


TFA 


2.75[nm] 








4.3171 


0.0986 


0.7077 


103.41 


NA 


2.75[nm] 


1 00 


0.3083 


3.4797 


0.0990 


0.9842 


104.93 


NA 


2.75[nm] 


1 000 


3.083 


2.3688 


0.1003 


1.8842 


111.71 


NA 


2.75[nm] 


5 000 


15.415 


1.7450 


0.1034 


3.1169 


127.69 


NA 


2.75[nm] 


10 000 


30.83 


1.5082 


0.1061 


3.8457 


141.13 


86.11 


2.75[nm] 


50 000 


154.15 


1.0202 


0.1179 


6.0508 


202.04 


163.9 


2.75[nm] 


100 000 


308.3 


0.8416 


0.1267 


7.2226 


248.1 


216.3 


2.75[nm] 


500 000 


1541.5 


0.5215 


0.1588 


10.501 


432.06 


411.80 


2.75[nm] 


1 000 000 


3083 


0.4226 


0.1784 


12.203 


559.92 


543.37 


2.75[nm] 


5 000 000 


15415 


0.2598 


0.2396 


17.070 


1044.7 


1034.4 



Table 9. Ground state chemical potential /i g , maximum value of the wave function 
<f> g (0, 0) and root mean square sizes r rms , z rras versus the interaction coefficient k 3 — k 
in 3d with cylindrical symmetry under case II. 
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compares the numerical ground state solution (i.e. numerical solution of (5.9)) and 
the Thomas- Fermi approximation (TFA) in (5.18). 

From Figures 4-5, and Tables 8-9, we can see that the chemical potential fi g and 
the root mean square sizes r rms , z Tms will increase when the interaction coefficient 
k 3 = k ( i.e. the number of atoms in the condensate) is increasing. On the other 
hand, the peak of the ground state solution 9 (O,O) will decrease. Figure 6 and 
Table 8-9 show that the Thomas-Fermi approximation are accurate for the chemical 
potential and ground state wave function near the origin only when k is very big, but 
gives poor approximation when k is intermediate or in the tail of the wave function. 

6.4 Application to compute excited states of GPE 

Suppose the eigenfunctions of the nonlinear eigenvalue problem (3.2) under the 
constraint (3.3) are 

±0 9 (x), ±01 (x), ±02 (x), 

whose energies satisfy 

E K (<f> g ) < E^fa) < E K (<j> 2 ) < ■ ■ ■ . 

Then <pj is called as the j-th excited state solution of the GPE (2.13). In fact, 
(f) g and 4>j U — 1)2, •••) are critical points of the energy functional E K ((f>) under 
the constraint (3.3). In Id, when V(x) — \ is chosen as the harmonic oscillator 
potential and K\ = 0, the excited states are given [33]: 

= (2>j!)~ 1/2 Q 1/4 Hj ( X ), 

H = M0j) = E o((f>j) = (j + £) , j = 1, 2, • • • ; 

where Hj(x) is the standard j-th Hermite function [33]. Here we show numerically 
that the algorithm (3.9) can also be applied to compute any j-th excited state of 
GPE with K\ > provided that we start with the above j-th excited state as initial 
data for k± > small and use a continuation technique for k± > bigger, i.e. use the 
numerical solution of the j-th excited state solution for a small K\ as initial guess for 
computing the j-th excited state of a larger K\. When the algorithm (3.9) is applied 
to compute j-th (j is an odd integer) excited state in Id, due to these functions are 
odd function, the finite element subspace in Section 4 should be replaced by 

V g h = {0 ft (r)GC([O, J R])|0 /l (r)| [rj , r3+l] GP 1 ([r,,r, +1 ]), < j < M - 1, 

0» = 4>\R) = o}. 

For simplicity, here we only report numerical results in Id for the first four excited 
states of GPE for different K\ > 0. Table 10 shows the energy E K (<f>j) and chemical 
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0.5000 


1. 


.5266 


3. 


.5966 


6.5527 


10.369 


19.070 


30.259 




1.5000 


2 


.3578 


4. 


3442 


7.2802 


11.089 


19.784 


30.971 




2.5000 


3 


.2590 


5, 


1479 


8.0432 


11.833 


20.512 


31.691 




3.500 


4. 


.1919 


5. 


.9901 


8.8349 


12.597 


21.252 


32.419 


/i 4 


4.500 


5. 


.1424 


6. 


.8598 


9.6501 


13.381 


22.005 


33.157 



Table 10. Energy and chemical potential of the ground state and first four excited 
states of GPE in Id (d = 1). 

potential fij (j = #,1,2,3,4) of the ground state and first four excited states of 
GPE in Id for different K\. Figure 7 plots the first four excited wave functions 4>j(x) 
(j = 1, 2, 3, 4) versus x for different k±. 

From the results in Table 10 and Figure 7, we can see that the algorithm (3.9) 
can be applied to compute the excited states of GPE (2.13). We observed from our 
numerical results in Table 10 that for any fixed Ki > 

E K (<P g ) < E^fa) < E K (<p 2 ) <■■■, n g < ii x < m < ■ ■ ■ . 

This implies that the eigenvalue of (3.2), fx g , corresponding to the minimizer of 
the energy functional E K ((f>), <p g , is the minimum chemical potential among all the 
eigenvalues of (3.2). Furthermore, we have 

lim |4M = 1, lim ^ = 1, j = 1,2,3,4. 

Kl ^+oo E K ((f) g ) Kl ^+oo y bg 

A rigorous mathematical justification of these numerical observations is under fur- 
ther study. 

7 Conclusions 

Ground state solution of time-independent Gross-Pitaevskii equation of Bose-Einstein 
condensation at zero or very low temperature is computed by directly minimizing 
the energy functional under a constraint through the finite element discretization. 
We begin with the 3d Gross-Pitaevskii equation, scale it to obtain a three-parameter 
model, show how to reduce it to 2d and Id GPEs. The ground state solution is for- 
mulated via minimizing the energy functional under a constraint. The finite element 
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approximation for Id, 2d with radial symmetry and 3d with spherical symmetry and 
cylindrical symmetry are presented in detail and approximate ground state solutions, 
which are used as initial guess in our practical numerical computation, are provided 
in two extreme regimes: very weak interactions and strong repulsive interactions. 
Numerical results are reported in Id, 2d with radial symmetry and 3d with spher- 
ical symmetry and cylindrical symmetry for condensation with repulsive/attractive 
interparticle interactions and atoms in it ranging up to millions to demonstrate the 
novel numerical method. Our numerical results show that the Thomas-Fermi ap- 
proximation are accurate for the chemical potential and ground state wave function 
near the origin only when k is very big, but gives poor approximation when k is in- 
termediate or in the tail of the wave function. Furthermore extension of our method 
to compute the excited states of GPE is also presented. In the future we plan to 
study physically more complex systems based on this ground state solution solver. 
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Figure 1: Ground-state condensate solution in Id (d = 1). a). Wave function 
<t> g {r) versus r for Kl = -12.5484, -6.2742, -2.5097, 0, 3.1371, 12.5484, 31.371, 
62.742, 156.855, 313.71, 627.42, 1254.8 (in order of increasing width), b). Chemical 
potential \x g versus Ka- 




b). K 2 
Figure 2: Ground-state condensate solution in 2d with radial symmetry, a). Wave 
function <f> g (r) versus r for « 2 = -5.8, -5.5, -4.5, -2.5097, 0, 12.5484, 62.742, 
313.71, 1254.8 (in order of increasing width), b). Chemical potential \i g versus K d 
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a). ' r b ). log(K 3 ) 

Figure 3: Ground-state condensate solution in 3d with spherical symmetry, a). 
Wave function <j> g {r) versus r for « 3 = -7, -6.2472, -3.1371, 0, 3.1371, 12.5484, 
31.371, 125.484, 627.42, 3137.1, 31371 (in order of increasing width), b). Chemical 
potential \i g versus n d . 
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Figure 4: Ground-state solution in 3d with cylindrical symmetry under case I. Con- 
densate wave function on two lines for k 3 = —4.705, —1.881, 0, 3.762, 9.405, 18.81, 
37.62, 94.05, 188.1, 376.2, 940.5, 1881, 3762 and 15048 (in order of increasing width): 
a). On the line z = {4> g (r, 0)); b). On the line r = {4> g {0, z)). Surface Plots of the 
condensate wave function (f) g (r,z): c). k 3 = 15048 (N = 800 000); d). k 3 = 188.1 
(N = 10 000). 



30 



0.1 0.2 0.3 0.4 0.5 

a). 



o 

b). 



10 



20 30 40 

z 




Figure 5: Ground-state solution in 3d with cylindrical symmetry under case II. 
Condensate wave function on two lines for k 3 = 0, 0.15415, 0.6166, 1.5415, 3.083, 
6.166, 15.415, 30.83, 61.66, 154.15, 308.3, 924.9, 1541, 3083 and 15415 (in order 
of increasing width): a). On the line z = ((f> g (r, 0)); b). On the line r = 
(<j) g (0,z)). Surface Plots of the condensate wave function <f> g (r,z): c). k 3 = 3083 
(N = 1 000 000); d). k 3 = 30.83 (N = 10 000). 
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Figure 6: Comparison between the numerical ground state solution and the Thomas- 
Femi approximation in 3d with cylindrical symmetry. ' — ': Numerical solution of 
(5.9), '+ H — h': Thomas-Fermi approximation (5.18). For case I: a). k 3 = 15048 
(N = 800 000); b). k 3 = 188.1 (N = 10 000). For case II: c). At the line z = 0, i.e. 
(j) g (r, 0); d). At the line r = 0, i.e. 0(0, z). 



32 




c). x d). x 



Figure 7: Excited states of GPE with harmonic oscillator potential in Id for K\ = 
0, 3.1371, 12.5484, 31.371, 62.742, 156.855, 313.71 (in order of increasing width), 
a). First excited state 0i(x) (odd function); b). Second excited state 02 (x) (even 
function); c). Third excited state 03 (x) (odd function); d). Firth excited state 4>±{x) 
(even function). 
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